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INTRODUCTION 

Randomly  irregular  waves  are  difficult  to  incorporate  into  engineering 
design.  The  response  of  the  structure  is  usually  quite  complicated  and  often 
involves  other  factors  than  wave  action.  Thus  it  is  frequently  not  possible 
to  analytically  determine  the  statistical  characterization  of  the  response 
directly  from  the  statistical  properties  of  the  waves.  In  problems  of  this 
type,  simulation  techniques  have  often  been  the  only  successful  method  for 
determining  solutions.  These  techniques  have  been  used  in  a  wide  range  of 
problems  in  physics,  operations  research,  and  other  fields,  wherever  random 
factors  were  involved  in  a  complicated  interaction  with  other  factors. 

Basically,  simulation  techniques  are  procedures  whereby  artificial  data 
having  imposed  statistical  properties  is  generated  by  sane  computational  means. 
Usually  this  is  done  in  a  digital  computer.  The  artificial  data  is  fed  into 
the  problem  and  the  response  calculated.  By  doing  this  with  enough  data,  the 
equivalent  of  many  years,  or  even  centuries,  of  experience  with  the  problem 
can  be  produced.  Such  factors  as  maximum  response,  or  the  number  of  times 
some  critical  value  is  attained,  can  be  determined  by  inspection  or  by  moni¬ 
toring  the  output  with  the  computer.  Simulation  techniques  have  the  advantage 
of  working  for  fairly  complicated  situations,  but  the  disadvantage  of  often 
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requiring  sizeable  amounts  of  computer  time. 

Simulation  procedures  have  not  been  used  extensively  in  coastal  engineering 
and  ocean  wave  problems  although  several  of  the  oil  companies  have  used  the 
techniques.  The  following  study  was  undertaken  to  make  the  procedures  more 
available  to  the  engineer  working  with  ocean  wave  problems  and  to  investigate 
possible  ways  to  increase  the  efficiency  and  the  realism  of  the  ocean  wave  and 
force  simulations  produced. 

The  conventions  used  in  defining  spectral  density  are  not  standardized. 
Differences  of  tt  and  2  show  up  in  various  papers  depending  on  whether  one¬ 
sided  or  two-sided  spectral  densities  are  used  and  whether  frequencies  are 
expressed  in  radians  or  cycles  per  unit  time.  All  derivations  in  the  following 
analysis  will  be  based  on  the  two-sided,  cycles-per-unit-time  spectral  density 
relations.  These  will  be  converted  to  one-sided  relations,  where  appropriate, 
by  multiplying  by  2  and  taking  the  integration  from  zero  to  infinity  instead 
of  from  minus  infinity  to  plus  infinity.  The  exact  mathematical  definitions 
are  given  in  the  table  of  notation  at  the  back  of  the  paper. 

The  sea  su.'face  elevations  will  be  assumed  to  be  a  stationary,  ergodic 
stochastic  process  produced  by  the  addition  of  many  infinitesimal  wavelets 

each  with  a  random  phase.  By  the  usual  random  theory  of  ocean  waves,  this 

2  3 

leads  to  a  Gaussian  process.  * 

2 

‘Pierson,  W.  J.,  Jr.,  "The  Representation  of  Ocean  Surface  Whves  by  a  Three- 
Dimensional  Stationary  Gaussian  Process,"  New  York  University,,  New  York,  195^» 
^Kinsman,  Blair,  Wind  Waves.  Prentice-Hall  Inc.,  Englewood  Cliffs,  N.  J., 

1965. 
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SIMULATION  OF  SEA  SURFACE  ELEVATIONS 


Two  basic  methods  for  simulating  ocean  wave  processes  were  studied.  These 

were  (a)  by  wave  superposition  and  (b)  by  linear  filters.  Each  method  has  its 

advantages  and  disadvantages.  Both  techniques  seek  to  produce  a  mean-zero, 

Gaussian  stochastic  process  which  has  an  initially  specified  function  as  its 

spectral  density.  This  initially  specified  function  will  be  called  the  target 

spectral  density,  while  the  spectral  density  estimated  from  the  simulated  time 

series  data  will  be  called  the  realized  spectral  density.  The  target  spectral 

4  5  6 

density  may  be  specified  by  a  theoretical  curve  *  *  or  by  the  discrete  tabu¬ 
lation  of  a  spectral  density  estimated  from  actual  ocean  wave  recordings. 

Simulation  by  Wave  Superposition.  -  Let  the  target  spectral  density  be 
denoted  by  s(f)  and  suppose  that  F  is  a  frequency  in  cycles  per  second  such 
that  s(f)  is  essentially  zero  if  f  is  greater  than  F.  Let 
0  s  f 0  <  f  1  <  f£  K  •••<%  =  ? 
be  a  partition  of  the  interval  (0,F)  and  define 

Af„  =  fn  -  fn-l  (1) 

K-  <fn+  W/2  (« 


^Kitaigorodskii,  S.  A.,  Application  of  the  theory  of  similarity  to  the 
analysis  of  wind  generated  wave  motion  as  a  stochastic  process,  Izv.  Akad. 

Nauk  SSSR  Ser.  Geofiz..  1,  105-117;  English  TransL,  1,  73-80,  1962. 

^Bre tschneider,  C.  L. ,  A  one  dimensional  gravity  wave  spectrum:  Ocean 
Wave  Spectra.  Prentice-Hall  Inc. ,  New  York,  1963. 

^Pierson,  W.  J.,  Jr.  and  L.  Moskowitz,  A  proposed  spectral  form  for  fully 
developed  wind  seas  based  on  the  similarity  theory  of  S.  A.  Kitaigorodskii: 
Jour.  Geophysical  Res.,  vol.  69,  no.  24,  pp.  5181-5190,  1964. 
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The  quantity  ^  (t)  given  by  the  formula 

*1  (t)  :  2  ZL  )J s(fnl  A  fn  cos  (knx  -  2  rr  f Rt  +  §>„>  .  (3) 

n=l 

where  <f>n,  n  z  1,2,  are  independent  random  variables  distributed 

uniformly  over  the  interval  (0,2  rr  )  and  kn  is  defined  by  the  relation 

(2  tt  fn)2  =  kng  tanh  knd  ,  (4) 

will  approximate  a  Gaussian  stochastic  process  with  zero  mean  and  spectral 

density  s(f).  The  symbol,  d,  denotes  water  depth  and  g  is  the  acceleration 

due  to  gravity.  A  standard  subroutine  is  available  on  most  computers  to 

generate  independent  uniform  random  numbers  that  may  be  used  for  <pn* 

anDroximation  improves  as  N  increases  and  max  A  fn  decreases. 

1<  nSN 

The  simulation  is  based  on  the  superposition  of  many  waves,  each  having  a 
random  phase  and  an  amplitude  consistent  with  the  energy  in  the  target  spectral 
density  at  that  frequency.  Fundamentally  it  is  just  the  finite-difference 
approximation  to  the  psuedo- integral  representation  for  ocean  waves. 

r°° 

n.  (t)  :  2  j  y/s(f)df  cos  (kx  -  2  rr  f t  +  $  )  (5) 

J0 

One  is  tempted  to  set  Afn  eqial  to  F/N  and  thus  use  an  equal  spaced 
subdivision  of  the  interval  (0,F).  However,  this  results  in  *\(t)  repeating 

A 

itself  exactly  with  period,  l/fj  .  There  are  several  ways  to  avoid  the 

^Pierson,  W.  J.,  Jr.,  "Wind  Generated  Gravity  Waves,"  Advances  in  Geophysics. 
vol,  2,  pp.  93-l?8>  Academic  Press,  New  York,  1955* 

®Brown,  L.  J. ,  "Methods  for  the  Analysis  of  Non- Stationary  Time  Series 
with  Applications  to  Oceanography,"  Hydraulic  Engr,  Lab.  Rep.  HEL  16-3. 

University  of  California,  Berkeley,  1967* 
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periodicity.  One  way  is  to  select  the  set  of  fn  values  with  a  random  number 
table.  Another  way  is  based  on  the  cumulative  spectrum,  defined  as 

S(f)  =  7./  s(f‘)df'  .  (6) 

The  quantity,  s(fn)  A  f  ,  is  approximately  the  same  as  [s(l‘n)  -  S(fn_^)]  / 2  . 
Hence  (3)  can  be  written  as 

N  _ 

q(t)  =  /T  X  /sTfn)  -  S(fn.i)  cos  (knx  -  2trfnt  +  $^)  (7) 


n=i 


The  periodicity  is  avoided  if  the  set  of  fn  values  are  chosen  to  make 

S(fn)  -  S(fn_^)  constant,  say  equal  to  a c  ,  for  all  n  values.  Then  (7)  becomes 


n(t)  s  /Ta  ^  cos  (knx  -  2Trfnt  +  $n) 
nrl 


(8) 


with  fn  defined  as  the  solution  of 

S(fn)  a  (n/N)  S (oo)  .  (9) 

This  corresponds  to  an  equal  subdivision  of  the  energy  coordinate  axis  for  the 
function  S(f)  .  Equation  (9)  is  particularly  easy  to  solve  if  the  Bretschneider- 
Pierson  spectral  density  is  used  as  a  theoretical  model.  This  model  has 
the  form 

-B/f4 


s(f)  =  ~  e" 
f5 


(10) 


and  is  directly  integrable 
S(f) 


:  2 

L  s5  2  L  Jo 


A  .-B/f" 
2 


(ID 


Hence  S(t»)  z  A/2  and  the  solution  of  (9)  is 


.1 


1/4 


log. 


(!)J 


(12) 
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The  constants  A  and  B  can  be  deduced  from  Pierson's  results  although  it  id 
necessary  to  be  careftil  about  the  distinctions  between  one-sided  and  two-sided 
spectral  densities  and  between  radian  and  cycles-per-second  frequencies.  The 
above  equations  are  based  on  a  two-sided  spectral  density  and  cycles-per-second 
frequency. 

The  function  S(f)  can  also  be  computed  from  measured  spectral  densities 
by  numerical  integration  and  graphical  determination  of  the  fn  values.  The 
list  of  fn  frequencies  is  then  read  into  the  simulation  computer  program  as 
input  data. 

The  simulation  for  sea  surfaces  having  a  directional  spectral  density, 
s(f, 6),  is  directly  analogous  to  (3)  and  is  just  the  discrete  analogue  to  the 
psuedo  integral 


n.(x*y*  t)  = 


oo 


/s(f,S)dfd6  cos(kx  cos©  +  ky  sin©  -2rrft-f-^>) 


(13) 


That  is 


N  M 


n.(x,y.t)  =  2  21  21  ^  8tfn»9m)  AfnA0m  cos(knx  cos  k^y  sin©^- 2  Trf^t+ 

n=l  mrl 

(14) 

where  the  interval  (0,2tt)  have  been  subdivided 
0  =  6q<  ©x<  ©2  <  •••  <  ©K  =  2TT 


and 


as  =  e  -em , 

m  m  m-i 


(15) 

e.  =  (»m+  e.-i  )/2  (i6) 

Equation  (14)  may  be  computed  for  more  than  one  space  location  using  the 
same  uniform  (0,2n)  independent  random  nunbers  for  each  computation.  If  N  and  M 

A 

are  large  enough,  the  simulated  values  of  ^(xj^yj^t)  ,  k  =  1,2,3»****  will 
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maintain  approximate  the  correct  intercorrelations  and  cross-spectral 
densities. 


Simulation  by  Linear  Filtering.  -  Let  Xp  x y  Xy  •••  ,  be  any  initial 

sequence  of  numbers  and  let  a_jj,  a  ^  p  a  *  *  *  ^  i»  aN  any  fixed 

sequence  of  constants.  Then  the  sequence  obtained  by 

N 

yk  =  ]>[  anxk-n  k  r  N  +  l,  N  +  2,  N  +  3,  •••  (17) 

nr  -  N 

is  called  the  output  obtained  by  applying  the  digital  filter 
lan  ,  n  s  0,  ±1,±2,  •**,±N }  to  the  initial  sequence  |xn,  n  s  1,  2,  •••]. 
The  basic  problem  (the  design  of  the  digital  filter)  is  the  determination  of 
the  values  of  aR  which  yield  yn  having  particular  desired  properties.  The 
procedure  for  designing  a  digital  filter  is  easy  to  follow  if  certain  basic 
relations  from  Fourier  transform  theory  are  kept  in  mind. 

Fundamental  Fourier  Relations.  -  A  linear  integral  operator  with  kemal, 
k(t),  acting  on  the  input,  x(t),  to  produce  an  output,  y(t),  m^y  be  written 

y(t)  =  fk(r)x(t-r)dr  .  (18) 


The  input  and  output  will  be  assumed  real;  hence,  the  kemal  function  is  also 


real.  The  Fourier  transform  of  the  kernal, 


K(f ) 


■/' 

-o© 


-i2nfr 

e  k(T)d*  , 


(19) 


is  called  the  system  function  of  the  linear  operator.  By  Euler's  relation, 
-i© 

e  =  cos  0  -  i  sin©  , 

SO  K(f)  can  be  also  written  in  terms  of  trignonetric  functions  as 

oo  ,  oo 

k(T)  cos  2nfrd*-  i  /  k(r)  sin  2n,ftd‘t 

-oo 


K(f ) 


■l 


oo 


(20) 


(21) 
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Relative  to  the  variable,  f,  the  real  part  of  K(f)  is,  this,  symmetric  about 


f  z  0  ,  and  the  imaginary  part  of  K(f)  is  stow  symmetric.  That  is 

<ft[K(f)]  r  #[K(-f)]  (22) 

j[K(f)]  r-Jl[K(-f)l  (23) 

The  two  basic  spectral  density  interrelationships  between  x(t)  and  the  y(t) 
o 

given  by  (18)  are 

:yy(f)  r  f  K(D|2  v(f)  (24) 

.^(f)  =  K(f)  .^(f)  (25) 


A  further  relation  is  that  y(t)  will  have  mean  tero  and  be  normally  dis¬ 
tributed  if  x(t)  has  these  same  properties. 

The  digital  filter  may  be  written  as  a  linear  integral  operator  if  Dirac 
delta  functions  are  introduced.  Suppose  that  and  x^  are  related  to  y(t) 
and  x(t)  by  the  following  equivalence. 

yk  =  y(kat)  ,  k  :  0,1,2, 3»*"  (26) 

Xn  =  x(nAt)  ,  n  s  0,1,2, 3, •••  (27) 

Then  the  digital  filter  in  (17)  may  be  rewritten 
N 

y(t)  r  ^  an  x(t  -  n  At) 
ns-H 


r°°N  I 

s  |  TV  an  £(f  -  n  M)1  x(t-  T)dT  ,  for  t  s  0,  At,26t,*** 

i<jo  nr- N 


Bendat,  J.  S.  and  Piersol,  A.  0., 

John  Wiley,  New  York,  1966,  pp. 96-99.  eqe.  (3.137)  end  (3.138). 
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(29) 


Thus  the  kemal  of  the  digital  filter,  which  will  be  denoted  by  k(T),  is 

N 

k(T)  r  y  »n  S(T  -  nAt) 
nr-N 

and  the  system  function  (derived  in  appendix  I)  is 
K(f )  r  f  e“i2nft  k(  t )  dr 

J-0O 


H  •  N 

:  ig  +  2  V  Ajj  cos  (ntrf/F)  -  i2  /  Bn  sin  (nTrf/F)  (30) 

n=-N  nr-N 


where 

Ajj  r  l/2  f  an  +  a_n] 

®n  -  1/2  £an  ~  a-nJ 

F  :l/2At 


(3D 

(32) 

(33) 


Equation  (30)  is  the  key  to  the  design  of  a  digital  filter  to  approximate 
an  arbitrary  linear  integral  operator  having  a  real  kernal.  Basically  one  has 

A 

only  to  make  K(f)  have  the  same  shape  as  the  system  function  K(f)  for  the 
arbitrary  operator.  That  is,  the  An  and  need  to  be  determined  so  that 

N 

(R[K(f)]  -  Aq  +  2  ^  Aj,  cos  (nTTf/F)  (34) 

nrl 


N 

j  [K(f)]  ~  2  2  Bn  sin  (n7Tf/F)  (35) 

nrl 

The  right  hand  side  of  these  equations  have  period  2F.  If  N  is  large  enough, 
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the  approximation  can  be  made  quite  good  between  -F  and  +  F  t  but  outside  that 

A 

interval  K(f)  will  be  repeated  with  period  2F.  However,  if  F  is  large  enough 
so  that  the  response  of  the  filter  to  frequencies  greater  than  F  is  of  no 
importance  in  the  intended  application,  then  the  periodicity  of  K(f)  will  cause 
no  difficulty.  Since  F  r  l/2At,  making  F  large  is  equivalent  to  moving  the 
digital  filter  along  a  sequence  with  tighter  spacing  on  the  time  axis. 

The  constants  An  and  Bn  may  be  determined  by  the  usual  procedures  for 
fitting  a  Fourier  series 


(36) 


(37) 


(Standard  subroutines  are  available  at  most  digital  computer  installations  to 

make  these  computations  quickly  and  easily.  The  new  fast  Fourier  methods  are 

10 

particularly  appropriate  here  .) 

Once  the  Ap  and  Bp  are  computed,  the  coefficients  an  follow  immediately 
from  (31)  and  (32) 


a0  '  A0 

an  -  ^n  +  ®n  (38) 

for  n  =  1,2,3»  ***»N 


^Cooley,  J.  W.,  and  Tukey,  J.  W.,  "An  Algorithm  for  the  Machine  Calculation 
of  Complex  Fourier  Series,"  Jour,  of  Math,  of  Computations.  April,  19&5* 
pp.  297-301. 
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Simulating  r^(t)  at  One  Soace  Location.  -  The  random  input,  x(t),  is 
called  white  noise  if  the  spectral  density  of  x(t)  is  unity.  White  noise  may 
be  approximated  by  a  sequence  of  independent  normally  distributed  random 
variables  x^^tX^, •  •*,  each  of  which  has  mean  zero  and  unit  variance.  Since 
subroutines  for  generating  independent  normal  random  variables  are  available 
for  most  cocputers,  white  noise  is  a  convenient  input  for  simulation  by  a  linear 
integral  operator  (18)  or  its  digital  approximator,  (17). 

Suppose  it  is  desired  to  produce  a  simulated  sea  surface  elevation,  (t) , 
which  has  spectral  density,  s(f).  This  can  be  produced  from  (18)  with  white 
noise  input  if 

K(f )  =  si  s(fj  .  (39) 

For  if  rj(t)  is  formally  identified  with  the  y(t)  produced  by  (18) ,  then 
(24)  gives 

sn(f)  =  |K(f)\2  sxx(f)  =  s(f)  (40) 

It  remains  only  to  determine  the  digital  filter  constants,  an ,  so  that  the 

A  I 

digital  filter  system  function  K(f)  closely  approximates  K(f)  =  V  s(f)  . 

This  is  achieved  by  the  Fourier  series  fitting  procedure  indicated  by  (34)  to 
(37).  Since  /s(f)  has  no  imaginary  part,  only  need  to  be 

determined.  The  cutoff  frequency,  F  ,  may  be  any  convenient  value  such  that 
s(f)  is  essentially  zero  for  higher  frequencies.  However  computations  are 
simplified  if  F  corresponds  to  a  At  interval  which  is  an  exact  fraction  or 
multiple  of  the  time  scale  unit. 

The  introduction  of  a  directional  spectral  density,  s(f, ®)»  as  the  target 
produces  no  real  complications  if  one  is  interested  only  in  one  space  location. 
All  that  is  necessary  is  to  integrate  out  ®  .  That  is, 
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3(f) 


(41) 


f27T 

s(f,  B  )  dd 

K 

and  then  proceed  as  in  (39)  and  following. 

Simulating  Several  Simultaneous  Time  Series.  -  The  basic  difficulty  with 
simultaneous  simulation,  is  that  the  individual  simulations  have  to  produce  the 
intercorrelations  or  interspectral  densities  between  the  various  series  which 
have  been  previously  specified.  One  procedure  which  maintains  these  inter¬ 
relationships  is  as  follows. 

Suppose  M  time  series  are  to  be  simulated.  These  will  be  denoted  by 
y1(t),  y2(t),  ••*,  y^(t).  The  simulations  will  be  developed  from  M  independent 
inputs,  x^(t),  x2(t),  •••,  Xy (t ) .  The  idealized  simulation  for  the  m-th  time 
series  will  be  given  by 


m 


^  =  Z 


oo 


knj(T)  Xj(t  -  “0  dT  ,  for  m  =  1,2,»**,M 


(42) 


>1  -oo 

(The  integral  operators  will  be  replaced  by  digital  filter  approximations  in 
the  actual  computations.) 

Let  smr(f)  represent  the  cross  spectral  density  between  ym(t)  and  yr(t). 

The  kernals  kmj(T),  or  the  corresponding  system  functions  ^(f),  can  be 
selected  to  produoe  the  required  cross  spectral  densities.  The  key  relationship 
utilized  in  the  determination  of  the  K-*(f)  is  (see  appendix  II) 


snr^ 


Z  Vr)  KrJ<«  V/O 


r  r  1,2, * • *,m 
m  r  1, 2,  •  • *,M 


(43) 


where  it  is  assumed  that  ra>r  .  (Note  s^_(f)  =  s__(f)  ).  The  over  bar  in 

m  mr 

the  previous  equations  indicates  complex  conjugation.  If  (43)  is  expanded  into 
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a  system  of  equations,  with  independent  white  noise  inputs,  it  becomes 


su(f)  ; 

S21(f)  =  K21(f)  Kll(f) 

*22^  =  lK21^l  +  |K22^I 

s21(f)  :  K^Cf)  Kn(f) 

s32(f)  =  K31(f)  K2l(f)  -J-  K32(f)  K22(f)  m 

s33(f)  =  |K31(f)|2  +  |K32(0|24  |K33(f)|2 

etc. 

This  system  of  equations  can  be  solved  seqientially 


*u<«  = 

'J*  n(f) 

Ka(f)  r 

S21(f)/^SU(f> 

K22(f)  = 

^s22(f)  -  |K21^f^|  ] 

K31(f)  = 

831(f)/  /sn(f) 

K32(f)  : 

[»32<f)  -  K31(f>  IC21(f)J  '  K22(f) 

£5 

ii 

[s33(0  -  |Vf)|2-|K32(f)|2jl/2 

etc. 
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The  design  of  the  simulation  equation  is  completed  now  by  determining 

A 

^(f)  as  in  (30),  (34)  and  (35)  to  fit  Kj^f)  as  determined  by  (46).  This 
in  turn  determines  the  digital  filter  coefficients,  needed  to  approximate 

the  kernal  associated  with  the  system  function  ^(f)  •  Then  the  f '  ial  simu¬ 
lation  equation  reduces  to 

m  N 

ya(k  At)  r  ZZanmj  xj,k-n  *  nsl»2»,**»M  ^ 

jsl  ns4l 

where  „  for  n  r  1,2, 3*4, •••  is  the  j-th  generated  sequence  of  independent, 

Jt  n 

mean  zero,  unit  variance,  normal  random  variables. 

If  the  input  x^(t)  is  not  white  noise  the  system  of  equations  represented 
by  (43)  can  still  be  solved  sequentially  with  only  slight  additional  compli¬ 
cations. 

Simulating  Several  Sea  Surface  Elevations.  -  Let  ^(t),  Ij^t) 

represent  the  sea  surface  elevations  at  time  t  at  M  space  locations.  If 
yn(t)  is  replaced  with  1m(t)  in  the  development  in  the  previous  section, 
then  (46)  gives  the  required  digital  filter.  The  spectral  densities  s^f), 
s22(f),  s^^(f),  etc.,  are  all  set  equal  to  the  target  spectral  density,  s(f), 
for  the  sea  surface  elevations,  which  is  the  same  for  all  locations.  The  cross- 
spectral  densities  have  a  slightly  different  formulation  depending  on  whether 
the  sea  surface  spectrum  is  unidirectional  (say  directed  along  the  x-axis)  or 
directional.  In  the  first  case,  suppose  the  m  locations  are  x^^*  *  •  The 

cross-spectral  density  between  ^(x^t)  *}(x  »t)  is^ 

■*''*'Brown,  L.  J. ,  and  Borgman,  L,  E.,  "Tables  of  the  Statistical  Distribution 
of  Ocean  Wave  Forces  and  Methods  for  the  Estimation  of  Cp  and  Cj<,"  Wave  Research 
Report  rlEL  9-7.  Hydraulic  Engineering  Laboratory,  University  of  California, 
Berkeley,  Calif.,  1966,  Appendix  C, 


14 


fwas mm***!*** 


"I 


s^f)  =  [  cos  k(*r  -  \)  -  1  sin  *0^  -  %)]  ( 

In  the  second  case,  let  s(f,  ©)  be  the  directional  spectral  density  for  the 

sea  surface  elevations  and  (x1,y1),  (x2,y2)*  (x^y^)*  •••*  be  th®  M 

space  locations.  Then 

r  2tt 


s(f )  : 


s(f,e)  d& 


s^f)  =  s(f,  ©)  cos  {k(xr-xm)  cos5  + k(yr-ym)  sin©}  d© 

'0 

'2rr 

-i  s(f,6)  sin{k(xr-x_)  cos©  +  k(y  -  ym)  sin©}  d© 

J .  (<*9) 

J  0 

Distortions  of  ^  (t)  to  Produce  Skewness.  -  Steep  waves  in  nature  have 
skewed  surface  elevations.  That  is,  the  crests  tend  to  be  higher  above  mean 
sea  level  than  the  troughs  are  below.  The  simulation  procedures  previously 
given  produce  Gaussian  or  normally-distributed  sea  surface  elevations  which 
are  not  skewed.  Obviously  high  wave  conditions  violate  a  number  ox  assumptions 
in  the  statistical  theory  of  waves,  most  importantly  perhaps  the  assumed  linear 
superposition.  However,  until  a  successful  statistical  theory  of  waves  of  finite 
height  is  developed,  the  engineer  is  faced  with  doing  as  well  as  he  can  with  the 
existing  theory. 

The  present  simulations  can  be  made  to  appear  a  little  more  like  the  real 
sea  if  skewness  is  introduced.  One  way  this  can  be  done  is  as  follows.  The 

A 

Gaussian  sea  surface  simulation,  ^(t),  is  developed  as  before.  Then  a  normal 

12 

to  gamma  transformation  which  maintains  a  zero  mean  and  introduces  the  required 


i4Campbell,  G.  A.,  "Probability  Curves  Showing  Poisson's  Exponential 
Summation,"  Bell  System  Technical  Journal,  vol.  2,  1923,  pp.  95-113* 
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skewness,  s^  , 

Qx  =  (\2  -  l)/3 

q2  =  (a3 .  7 n.)/36 

Q3  r  (-31a4  -  ?^l2+  16)/810 

Q4  =  (9n5  +  256 n.3  -  433^/38880  (  50) 

Q5  :  (12H.6  -  243  rf  -  923  »T.2  +  14?2)/204120 

Q6  =  (-3753  ^ 7  -  4353  h,5  +  28951?  i3  +  289717  ?l)/l46966400 

6 

1  +  Z  «n  <«k  /  »” 

nrl 


is  made  to  introduce  more  peaked  crests  and  flatter  troughs.  Skewness  may  be 
introduced  in  other  ways.  The  above  is  just  one  convenient  and  thoroughly 
investigated  procedure.  The  whole  topic  deserves  considerable  further  research. 

How  much  skewness  should  be  introduced?  If  an  extensive  piece  of  actual 
wave  record  is  available  for  the  sea  surface  elevations  that  are  to  be  simu¬ 
lated,  one  scheme  of  calculation  would  be:  (a)  Compute  the  skewness  for  the 
wave  record 


sk 


3(t)  dt  / 


(5D 


(where  R  is  the  length  of  record),  (b)  Make  the  gamma-to-normal  transformation^3 
which  produces  an  unskewed  version  of  the  wave  record  with  the  same  zero  mean 
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Riordan,  John, 


"Inversion  Formulas  in  Normal  Variable  Mapping,"  Annals 


of  Mathematical  Statistics,  vol.  20,  1949,  pp.  417-425. 


xx  :  (-  ni  -+- 1)  /3 

*2  =  (?1s  -  'Is5  /J6 

x,  r  (-219 14  +  14H?  + 13)  /1620 

J  S3 

x4  =  (3993nf  -  152a  1+  119  n)  /3888O  (52) 

X-  =  (-67227  q6+  1707  ^  -  2041  *[ 2  -  3095)  /816480 

J  S  s  8 

X,  :  (10059417 1 J  -  179223  ll5  +•  271427  ^  3  +  215827  Y[  )  /146966400 

o  S  S  S  S 

6 

a<t)  r  1^  +2  X„  (sk/2)n 
n=l 

This  system  of  equations  gives  two  more  ten. «  than  were  contained  in  Riordan's 
paper  and  also  corrects  a  misprint  in  the  X4  formula.  Both  (50)  and  (52)  can 
be  truncated  to  fewer  terms  if  the  accuracy  requirements  permit,  (c)  Compute 
the  spectral  density  on  the  unskewed  version  and  make  a  normal  simulation, 

A 

r^(t).  (d)  Re-introduce  the  skewness  calculated  in  (51)  by  making  the  normal- 

to-gauma  transformation  given  by  (50). 

The  above  schene  needs  to  be  evaluated  against  actual  experience.  It 

represents  only  one  possible  procedure.  Undoubtedly  further  study  will  suggest 

others.  Perhaps  the  mode  of  the  spectral  density,  the  rms  wave  amplitude, 

14 

and  the  water  depth  could  be  used  to  get  Yl  /H  from  existing  graphs.  This, 

c 

in  turn,  might  provide  an  estimate  of  the  skewness  that  needs  to  be  introduced. 
14 

Isan,  R.  G. ,  "Stream  Function  Wave  Theory;  Validity  and  Application," 
Coastal  Engineering  Santa  Barbara  Specialty  Conference.  A5CE,  October,  1965, 

p.  282. 
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SIMULATION  OF  VELOCITIES  AND  ACCELERATIONS 


The  horizontal  components  of  the  water  velocity  and  acceleration  at  some 
fix*d  elevation  above  the  sea  floor  can  be  simulated  with  filters  using  either 
white  noise  or  the  sea  surface  elevations  as  basic  input.  A  unidirectional 
spectrum  and  a  white  noise  input  requires  the  same  procedures  previously 
developed  except  that  the  spectral  densities  for  the  horizontal  velocity  and 
acceleration  need  to  be  introduced. 


£>w(f) 


saa<f> 


(2TTf)2 


2 

cosh  kz 
sinh^kd 


(2nf)4 


cosh^kz 

sinh2kd 


3  (f) 

*1*1 


3  „n(f) 

’ll 


(53) 

(54) 


The  cross-spectral  density  between  the  horizontal  velocity  and  acceleration  is 
zero.  Hence  the  two  digital  filters  can  be  designed  separately  and  applied  to 
two  independent  white  noise  inputs. 

The  situation  for  a  directional  spectral  density  is  considerably  more 
complicated.  For  one  thing,  there  are  two  horizontal  components  for  both  the 
velocity  and  the  acceleration.  The  cross-spectral  densities  between  these 
components  are  not  all  zero.  Thus  the  procedure  for  simultaneous  simulations 
of  time  series  needs  to  be  introduced  (see  (42)  and  (43)).  The  required  cross- 

15  3L 

spectral  densities  are  listed  by  Wave  Research  Report  9-12. 

The  generation  of  velocities  and  accelerations  using  the  sea  surface 
elevation  as  input  is  based  on  (25)  and  the  cross-spectral  densities  of  velocity 


l^a 

Borgman,  L.  E. ,  "Tables  of  Ocean  Vote  Cross -Spectral  Formulas,"  Wave 
Research  Report  HEL  9-12.  Hydraulic  Engineering  Laboratory,  University  of 
California,  Berkeley,  Calif.,  196?. 
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and  acceleration  with  r^(t).  The  required  spectral  densities  are,  for  the 
unidirectional  spectrum 


1(x,t),  v(x+h,t) 


s  (f) 

Wx.t),  a(x+h,t) 


(f)  r  [(2TTf)  (cos  kh  -  i  sin  kh)]  s„a(f)  (55) 


[(2TTf) 


2  cosh  kz 
sinh  kd 


(sin  kh  +  i  cos  kh)]  s^(f)  (56) 


Hence  by  (25) 

K^f;*)  =  (27T  f )  (cos  kh  -  i  sin  kh) 


K  (f;z)  z  (27Tf)2  ^2§hjcz  (gin  kh  _  i  cog  kh) 
a  sinh  kd 


(57) 


(58) 


Thus  for  a  unidirectional  spectral  density,  the  system  functions  Ky(f;z)  and 
Ka(f;a)  do  not  depend  on  the  particular  sea  surface  spectral  density  that  is 
to  be  used.  The  digital  filters  can  be  determined  by  (3^)  and  (35)  without 


references  to  s^^(f)  .  This  yields  an  estimator  equation  that  may  be  written 


N 

v(k  At)  =  ^  a^  l|(k  A  t  -  n  A  t) 


nr-N 

N 


a(k  At)  z  aan*7.(k  At  -  n  At) 


(59) 


(60) 


ns-N 


Somewhat  similar  procedures  were  used  by  Reid  in  a  study  of  wave  forces^^5. 

The  corresponding  procedure  for  a  directional  spectral  density  is  analogous. 

15a 

The  cross-spectral  densities  may  be  obtained  from  Tech  Report  HEL  9-12.  The 
system  functions  intrinsically  depend  on  the  directional  spectral  density  for 
the  sea  surface. 

Both  the  unidirectional  and  the  directional  spectral  density  cases  can  be 

"^Reid,  R.  0.,  "Correlation  of  Water  Level  Variations  with  Wave  Forces  on 
a  Vertical  Pit  for  Nonperiodic  Waves,  Proc.  Sixth  Conf,  on  Coastal  Eng.,  pp.  7**9- 
786,  The  Engineering  Foundation,  Univ.  of  Calif. ,  1958. 
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handled  -airly  directly  with  the  procedure  for  simulation  by  wave  superposition. 
The  linear  wave  theory  formula  for  the  quantity  to  be  simulated  is  written  down 
using  as  wave  amplitude,  Jks($n)  A  fn  in  the  unidirectional  case,  and 
V 4s (?n,  0^  Afn  A0m  in  the  directional  case.  For  both  oases  a  random  phase 
is  inserted.  In  the  directional  case,  kx  in  the  linear  wave  theory  is 
replaced  with  kx  cos  8  +  ky  sin  6  and  the  whole  term  is  multiplied  by  cos  & 
for  the  x  component  and  sin  S  for  the  y  component.  Then  the  expression 

A  A 

is  sunned  over  all  fn  (and  9m  in  the  2-D  situation).  Examples  of  this  pro¬ 
cedure  are 

linear  theory ; 

v  =  a(2rrf)  cos  (kx  -  2* ft)  (6l) 

x  sinh  kd 


simulation  for  unidirectional  spectrum; 
N 


vx(t)  s^"  fin9nWn  (2tt fn)  C?f.h.  *n-  cos (knx  -  2  n  f nt -<•  $n)  (62) 

■  A  A  nil  b  A 


nal 


sinh  knd 


simulation  for  directional  spectrum:  v  (t)  z 

x 


N  M 

Z  v/4s(fn»5m)AfnAgm  (2^n)T""L  cos  0mcos(knx  cos^+k^  sin©m-2nfnt+§mn) 
nsl  m'1  (63) 


sinh  knd 


The  simulation  by  wave  superposition  has  the  one  grave  disadvantage  of  being 
time  consuming.  However  this  is  balanced  by  the  fact  that  simultaneous  simu¬ 
lations  of  several  quantities  automatically  maintain  the  proper  intercorrelation. 


SIMULATION  OF  WAVE  FORCES 


The  usual  formula  for  wave  forces  on  a  one-foot  section  of  vertical  piling 
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•* 


a 


(64) 


f  Cd|Dv|v|+Cm^  s£ 

can  be  applied  to  the  velocity  and  acceleration  simulations  developed  in  the 
previous  section  to  give  force  simulations.  This  procedure  could  be  followed 
to  simulate  the  total  force  on  a  structure  composed  of  vertical  cylinders. 
Simulations  could  be  made  for  the  forces  at  a  series  of  locations  down  each 
pile*  and  the  forces  could  be  numerically  integrated  to  give  total  force. 

Modifications  to  Give  Better  Agreement  with  Nature.  -  The  statistical 
theory  for  waves  is  based  on  linear  wave  theory.  But  linear  waves  have 
infinitesimal  amplitudes  which  do  not  extend  measurably  above  still  water 
level.  How  does  one  then  compute  velocities  and  accelerations  for  elevations 
above  this  level?  To  what  elevation  up  the  pile  does  the  integration  proceed 
if  it  is  not  carried  just  to  still  water  level?  If  the  velocities  and  accel¬ 
erations  are  simulated  from  the  sea  surface  elevations  at  the  pile,  the  second 
question  is  answered  immediately.  The  integration  is  carried  to  the  sea  sur¬ 
face.  The  firet  question,  however,  can  only  be  answered  with  approximations. 

A  conservative  approximation  for  velocities  and  accelerations  above  mean  sea 
level  is  to  use  the  same  formulas  that  bold  below  sea  level.  The  forces  are 
computed  for  z>d  and  ignored  if  (t)  <  z  .  An  alternative  approximation 
is  to  "stretch"  the  still  water  level  in  the  formula  up  to  the  sea  surface. 

In  this  procedure,  the  force  is  simulated  at  a  series  of  elevations  between 
the  sea  floor  and  still  water  level.  The  force  at  elevation  z  is  assumed 
to  act  actually  at  elevation  z[d+i^(t)j/d  .  That  is  the  force  displaced  upward 
or  downward  depending  on  whether  (t )  is  positive  or  negative.  This  pro¬ 

cedure  has  the  advantage  that  the  integration  for  total  force  can  be  carried 
to  z  =  d  in  terms  of  the  force  simulation  and  the  total  force  then  "stretched" 
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to  correct  for  i^(t)  . 

It  is  difficult  to  evaluate  the  accuracy  of  these  approxinations.  They 
both  obviously  fail  if  the  waves  get  very  high.  Yet  they,  at  least,  provide 
one  approach  to  the  problem.  The  ultimate  answer  is,  of  course,  a  statistical 
theory  for  waves  of  finite  height.  Until  that  is  available,  irregular  waves 
of  non-permanent  form  having  a  statistical  nature  can  only  be  treated  with  a 
"patched-up"  linear  statistical  theory. 

The  force  simulation  with  (64)  being  applied  to  (59)  and  (60)  and  then 
numerically  integrated  has  one  grave  disadvantage  in  computer  computations. 

It  is  time  consuming.  An  approximate  procedure  which  is  faster  would  be  better. 
After  all,  (64)  is  known  to  be  only  an  imperfect  approximation  of  the  actual 
force. ^  So  some  degree  of  further  approximation  is  not  particularly  objec¬ 
tionable. 

Approximations  to  the  Force  Formula.  -  A  faster  simulation  v^uld  result 
if  the  total  force  on  a  vertical  pile  could  be  expressed  as  a  single  operator 
acting  on  the  sea  surface  elevations  or,  alternatively,  acting  on  white  noise. 
The  inertial  force,  the  second  term  on  the  right  hand  3ide  of  (64),  can  be  so 
expressed.  If  a(t;z)  is  the  acceleration  at  elevation  z  on  the  pile  at 
time  t  and  k  (r  ;z)  is  the  kemal  for  a(t;z)  then 

cl 


a(t;z) 


*[(t  -  T  )  dr 


(65) 


Wiegel,  R.  L.,  Beebe,  K.  E.,  and  Moon,  James,  "Ocean  Wave  Forces  on 
Circular  Cylindrical  Piles."  Journal  of  the  Hydraulics  Division,  ASCE,  vol.  83, 
no.  HY2,  Proc.  Paper  1199,  April,  1957. 
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The  total  inertial  foroe  becomes  then 


total  inertial  force 


[(^(D+dj/d]  f  <*i! 

J  n 


OO 


ttD 

4 


ka(T  i*)  ^(t-T)dr  dt 


-co 


r°°  2  r d  1 

=  [(  n  (t)  +  d)/d]  ka(r;z)dt  ^(t-T)d* 

J-ooL  ?  J0  (66) 


This  eliminates  the  vertical  integration.  The  quantity  in  the  bracket  is  the 
new  total  kemal  that  must  be  approximated  with  a  digital  filter.  The  "stretching" 
correction  for  forces  above  still  water  level  has  been  introduced  into  the 
equation. 

The  drag  force,  the  first  term  on  the  right  side  of  (64),  does  not  simplify 
this  way.  If 

r  00 

v(t;z)  =  J  kv(T  ;z)  -  T)  dr  (6?) 

“  Oo 


is  substituted  into  the  expression  for  drag  force,  one  gets  for  the  total  drag 
force  on  a  vertical  pile 


fd  oa 

f°° 

[(  rj(t)  +  d)/d]  Cd^D  kv(r  ;z)  ^(t-t)dr 

1  kv(r  ;z)  ^(t-T)dx 

J0  J-oo 

■'-oo 

dz 


(68) 


The  absolute  value  interfers  with  the  shuffling  of  integrals  which  simplified  (fc£>). 

1?  ,  , 

In  an  earlier  paper,  the  term  v|v)  in  the  drag  force  formula  was  found 
to  be  approximated  very  nicely, at  least  relative  to  predicting  the  force  spectral 


Borgnan,  L.  E.,  "Spectral  Analysis  of  Ocean  Wave  Forces  on  Piling," 
Journal  of  Waterways  and  Harbors  Division.  ASCE,  vol.  93*  no.  WW2,  Proc.  Paper 
5247,  May  196?,  pp.  12^-156.  (See  p.151,  eq.(?4)). 
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density,  by 


/vj 

VI VI  =  V 


rms 


\j  (8/ tt) 


(69) 


Here,  v^g  is  the  standard  deviation,  or  root-mean- square  value  for  the 
velocity.  In  essence  |v|  has  been  replaced  by  v^g  \J (,8/tT)  .  The  effect 
of  the  approximation  is  easier  to  comprehend  if  the  linear  theory  versions  are 
introduced.  From  the  usual  formulas  for  velocity 

,2 


v|v|  =  fas  ssai-hil  cos  8  I  cos  8 1 
11  [  T  sinh  kd  J  1 


rms 


/Tv  :  ££SLMl2  HE  COS  0 

V  TT  (.  T  sinh  kd  J  y  2  TT 


(70) 


(71) 


Since  ^  Q/ztt  is  not  too  much  larger  than  unity,  the  approximation  replaces 
v | v{  with  a  cosine  oscillation  whose  maximum  amplitude  (at  9  =  0)  is  only 
slightly  greater  than  v|vj  .  Both  expressions  are  zero  at  6  s  TT/z  . 

The  analytic  basis  for  (69)  becomes  clear  if  one  asks  for  the  most  accurate 
linear  estimate  of  v|v| ,  assuming  v  is  normally  distributed  with  mean  zero 
and  standard  deviation  C  ~  vms  .  That  is,  what  is  the  constant  c  which 


minimizes 


Q  = 


oo 


[v|v|  -  cv] 


2  b-v2/2CT2 


-  Oo 


fzn  V 


dv 


(72) 


By  calculus,  Q  has  an  extremum  if 

C 


0  =  -  -2 
a  c 


r  00 
2/ ..i  e 


e-v2/ 2C72 

v^l  v|  - —  dv  -  c 


—  oo 


fzn  cr 


r.  dv] 


oo 


fzn  <r 


(73) 


or 


,.-v2/20-2 


r  *** 
fv2| 

J-  oo 


fzn  <r 


dv 


c  = 


cr /(8/tt) 


(74) 
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Further  analysis  shows  that  this  value  minimizes  Q  .  The  result  agrees  exactly 
with  (69) ,  hence  cv  is  the  best  linear  estimator  for  v|vj  in  the  least- 
square  sense.  The  procedure  can  be  extended  to  the  cubic  and  quintic  approxi¬ 
mations.  There,  one  seeks  the  values  of  the  constants  which  minimize 

2  .2/„^2 


<3  - 


o o 


-  Oo 


r  3 1  e“v  h 

V  V  -  C-.V  -  c,v^  fi - dv 

L  3  J  >T2rF  O' 


(75) 


or 


00 


[vM  - 


clv  "  C3V  "  C5V" 


-  00 


512  <r£l?ZL 
'Jz  tt  <r 


dv 


(76) 


The  solutions  which  minimize  (?5)  and  (76)  are  determined  by  differentiation 
with  respect  to  the  unknown  constants  to  be 


clv  +  C3V*  z  vrms  |/  (2 Itt)  v  +  (v3/vrms3  (?7) 

and 

Clv  +  c3v34-  03V5  =  |/(2/tt)  [(3vms/4)  V  +  (1/2  vrag)  v3-  (l/60  v^g)  v5] 

(7o; 

A  comparison  of  the  linear,  cubic,  and  quintic  approximations  with  v|v|  is 
shown  in  Fig.  1.  The  approximations  nave  been  :sade  dimensionless  by  changing 
to  the  variable 

X  =  (vhrm.)  (79) 

Then  the  three  approximations  become  in  terms  of  x 

(v!v^vras)  =  XM  ~  /^/rr)  x  (8°) 

V (2/T T)  [x  +  (1/3)  x3]  (81) 

[  (3/4)x  +(l/2)x3-  (1/60)*5]  (82) 
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Fig.  1.  Polynomial  approximation  to 
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An  examination  of  Fig.  1  shows  that  the  linear  approximation  does  a  fair  job 
for  j  x  |  <  2  ;  that  is,  if  the  velocity  is  within  plus  or  minus  two  standard 
deviations.  Since  according  to  normal  theory  this  happens  9 5i>  of  the  time, 
it  is  entirely  reasonable  that  the  linear  approximation  would  do  a  good  job 
on  estimating  a  curve  like  the  spectral  density  which  depends  on  the  average 
behavior  of  the  velocity  values. 

The  cubic  approximation  is  quite  accurate  out  to  more  than  three  standard 
deviations  and  the  quintic  is  good  to  almost  four  standard  deviations.  The 
quintic  approximation  is  only  slightly  better  than  the  cubic.  The  cubic 
approximation,  however,  is  quite  an  improvement  over  the  linear  expression. 

All  of  the  above  approximations  were  derived  under  the  assumption  that 
the  mean  value  of  v  was  zero  (i.e.,  there  was  no  steaty  current  present). 

The  same  analysis  can  be  carried  through  if  the  mean  of  v  is  not  necessarily 
zero,  although  the  mathematics  is  more  complicated.  Even  powers  of  v  must 
be  included  and  the  normal  density  has  to  include  the  parameter  for  the  mean. 

As  an  exanple,  the  least-square  linear  estimate  for  a  normal  v  with  mean  m 

2 

and  variance  O'  is  determined  by  minimizing 


e-(v-m)2/2CT2 

V  2  7T  <T 


dv 


(83) 


Let 


(84) 


e-t2/2  dt 

{ZTT 


The  least  square  solution  for  (83)  is 


(85) 
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+  V  .  «r2-»2)[2$(^-)  -lj  -2mO-<t>(i)  + 


+  {an[2$(i.)  -  l]  -*a  +(£))  v 


(86) 


Total  Force  on  a  Single  Vertical  Pile.  -  The  filter  concept  in  (18)  can 
be  extended  to  higher  orders 
oo 


y(t) 


oo  Q0 

I 

(t  )  x(t-  r)dr  + 

-  Oo  -  oo  Oo 

°°  r°°  r  OO 


*1 


r2)  xft-r^  x(t-r2)  drx  d  r2 


moo 

k3^Tl»  r2»  T3^  x(t"Tl)  x^t-  ^  x(t_7V  d7<idr2dr3 


-OO  -  OO  —  Oo 


(87) 


A  second-order  filter  would  consist  of  the  first  two  terms,  a  third-order  filter 
of  the  first  three  terms,  etc.  If  the  first-order  approximation  to  the  drag 
force  is  used,  the  difficulty  in  (68)  is  avoided  completely.  The  total  drag 
force  becomes  (approximately) 


Thus  the  quantity  in  the  bracket  is  the  kemal  for  the  estimation  of  total 
force.  The  system  function  for  the  total  drag  force  will  be  the  Fourier 
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transform  of  the  bracketed  quantity  or 


K 


(f) 

total  drag 


Similarly  from  (66) 

K  (f)  = 

total  inertial 


K  (f;z)  dz 
& 


(89) 


(90) 


and 


K  (f)  =  K  (f)  +  K  (f) 

total  force  total  drag  total  inertial 

The  system  functions  within  (89)  and  (90)  are  given  by  (57)  and  (58).  The 
integrations  in  (89)  ani  (90)  are  easy  to  make  and  the  result  is 


K  (f)  Z  ~~  Dv 

total  force  D  2S  rms 


(cos  kh  -  i  sin  kh)  + 


+  O,  JL  212 1  fart)2  (8in  kh  -  i  cos  kh) 
n  g  4  k 


(9D 


(92) 


This  system  function  can  be  approximated  with  a  Fourier  series  to  obtain 

A 

K  (f)  and  with  it  the  digital  filter  which  acts  on  the  sea  surface  to 

total  force 

simulate  the  total  force  on  a  single  vertical  pile. 

The  above  development  was  based  on  a  unidirectional  wave  spectrum  and  a 
linear  approximation  to  drag  force.  The  cubic  approximation  would  lead  to  a 
third-order  filter.  The  analysis  would  be  parallel  to  the  above  but  somewhat 
more  complicated. 

The  complications  cam  be  avoided  by  a  slightly  different  approach  which 
involves,  however,  many  more  Fourier  fits  to  determine  digital  filters.  Suppose 
the  digital  filters  for  velocity  and  acceleration  are  computed  at  a  series  of 
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points  up  the  pile  between  the  sea  floor  and  still  water  level. 

Then  the  force  at  the  k-th  position  is  (using  the  cubic  approximation), 
N 

fk(m&t)  z  a£^  f[(mAt  -  nAt) 

n»-N 


N  N  N 


IV  .  . 

Z  Z  Z 

n  z-N  n  =-N  n  =-N 


(mAt  -n*At) 


(93) 


with 


(k) 
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fn 
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:  t  f  nil  (v)  +  0,  JL  a(k) 

D  2g  V  7T  11113  m  g  4  an 


g 

(k)  (k)  (k) 

(k)  „  «  ^  Ho  %n'  &vn"  avn"' 


,(k)  C  JL  D  jCT 

f}n',n",nw  D  2g  V  9rr 


(94) 


(95) 


rms 


(k)  (k) 

and  a  and  a  are  the  digital  filter  coefficients  at  the  elevation  for 
vn  an 

the  k-th  position  on  the  pile.  The  total  force  may  be  computed  by  numerical 
integration  of  (93)  to  provide  a  cubic  digital  filter  for  total  force. 

N 

\r—  (x) 
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(97) 
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(k)  (k)  (k) 


af;n\n",n'"  =  ^  afn'  afn"  afn"'  A“k 


(98) 
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where  A  z^  is  the  appropriate  vertical  increment  for  the  k-th  position  on 
the  pile. 

The  total  force  can  be  "stretched"  to  correct  for  sea  surface  elevation 

variations  by  multiplying  F  (mAt)  by  Ln(m At)  +  dl /d  . 

total  L 

Total  Force  on  a  Chroup  of  Piles.  -  One  procedure  for  simulating  the  total 
force  on  a  group  of  vertical  piles  would  be  as  follows.  The  sea  surface 
elevations  are  first  simulated  for  some  central  position  within  the  group  of 
piles.  This  is  used  as  a  conmon  input  for  the  digital  filter  simulation  of 
the  force  on  each  of  the  piles.  That  is,  each  pile  would  have  a  different  h 
value  in  (57),  (58) *  and  (92),  which  would  be  the  x-direction  spacing  between 
the  central  position  and  the  pile.  If  desired,  the  sea  surface  elevation  at 
each  pile  could  also  be  simulated  from  the  sea  surface  elevations  at  the  central 
position.  All  these  together  could  then  be  used  to  determine  the  time  series 
for  total  force.  Of  course,  if  the  total  force  on  each  pile  is  not  going  to 
be  corrected  for  soa  surface  elevations  at  that  pile,  the  sea  surface  simu¬ 
lations  can  be  dispensed  with.  In  this  case,  the  force  digital  filter  coeffi¬ 
cients  for  each  pile  can  be  added  to  give  a  single  digital  filter  which  acts 
on  the  sea  surface  elevations  at  the  central  position  to  provide  as  output  the 
simulation  for  total  force  on  the  group  of  piles. 


AN  APPLICATION 


16 

Wave  forces  were  measured  near  Davenport,  California  in  1953.  The  values 
of  CD  and 

The  sea  surface  elevations  were  scanned  and  "well-defined"  waves  were  selected. 
The  wave  force  was  read  from  the  record  at  the  crest  and  at  the  still  water 
level  crossing  of  each  wave  profile.  The  height  and  period  for  each  wave  was 


C^  were  determined  for  a  nunber  of  waves  by  the  following  procedure. 
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used,  with  linear  wave  theory,  to  compute  the  water  velocity  and  acceleration 
in  the  vicinity  of  the  force  sensing  instrument  at  the  crest  and  zero  crossing 
phases  of  the  wave.  These  theoretical  velocities  and  accelerations  were  used 
with  the  measured  forces  to  determine  and  .  The  procedure  is  consid¬ 
erably  simplified  by  the  theoretical  relation  that  the  inertial  force  is  zero 
beneath  the  crest  and  the  drag  force  is  zero  beneath  the  zero-crossing. 

The  values  of  and  so  computed  showed  a  very  large  scatter.  There 
are  a  number  of  possible  reasons  why  this  might  happen.  The  forces  might  be 
inherently  random  so  that  CQ  and  Cyj  actually  do  vary.  Vibrations  or  other 
recording  difficulties  may  have  introduced  inaccuracies.  Perhaps  the  irregu¬ 
larities  in  the  wave  profiles  were  sufficient  to  invalidate  the  theoretical 
computations.  And  there  are  other  possibilities. 

As  a  means  of  estimating  the  effect  of  profile  irregularity  on  the  Cp 
and  Cjj  computations,  a  time  history  of  sea  surface  elevation  and  water  velocity 
and  acceleration  was  simulated  by  the  method  of  (7)  and  (62).  The  spectral 
density  for  the  recorded  waves  on  roll  10  of  the  Davenport  data  was  used  in  the 
simulation.  The  agreement  of  the  spectral  density  of  the  simulated  sea  surface 
and  the  Davenport  spectral  density  is  shown  in  Fig.  2.  The  velocities  and  accel¬ 
erations  were  simulated  at  the  same  elevation  above  the  sea  floor  as  was  used  in 
the  roll  10  measurements.  The  force  at  the  instrument  level  was  computed  from 
the  simulated  velocities  and  accelerations  using  (64)  and  a  Cp  value  of  1.0 
and  of  2.0  .  Hence  the  simulations  provided  a  sea  surface  time  history  and 
the  force  record  which  would  result  if  the  usual  force  formula,  (64),  were  to 
hold  exactly  with  the  specified  values  of  and  C^. 

Now  the  same  analysis  procedure  used  for  the  Davenport  data  was  applied 
to  this  sea  surface  and  force  record.  Well  defined  waves  were  selected  from 
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Roll  10  Actual 
Roll  10  Stimulated 

d  =  49  Feet 


0  0.04  0.08  0.12  0.16  0.20  0.24  0.28 

FREQUENCY  ( cycles /sec.) 

Fig.  2.  Simulated  and  actual  spectral  densities  for  Davenport  data,  roll  10. 


the  sea  surface  elevations.  The  forces  were  read  from  the  simulation  beneath 
each  wave  crest  and  wave  zero-crossing.  Linear  theory  was  used  to  compute  the 
velocity  beneath  the  crest  and  the  acceleration  beneath  the  profile  zero-crossing. 
These  values  were  combined  to  estimate  and  C^.  If  the  method  were  accurate, 
CD  z  1.0  and  z  2.0  would  be  the  result. 

The  scatter  diagram  for  the  computed  and 
There  is  considerable  scatter.  Fig.  4  provides  a  plot  of  Cq  ver3U3  Reynolds 
nunber  and  Fig.  5  gives  a  ranked  plot  of  Cjj  on  normal  probability  paper.  The 
figures  are  strikingly  similar  to  the  corresponding  results  for  the  Davenport 
data  and  suggests  that  the  wave  profile  irregularities  contributed  considerably 
to  the  CD  and  scatter  in  the  Davenport  results. 

CONCLUSIONS 

3.  Techniques  for  simulating  ocean  waves  are  most  applicable  when  the 
response  of  the  structure  is  conplicated  and  perhaps  involves  other  random 
environmental  factors  that  may  be  introduced  by  concurrent  simulations. 

2.  The  accuracy  of  the  wave  simulation  is  greatest  for  low  amplitude 
waves  and  decreases  for  large  steep  waves.  The  degree  of  loss  in  accuracy 
for  the  higher  and  steeper  waves  deserves  further  research. 

3.  Simulation  techniques  have  the  disadvarxage  of  being  time  consuming 
and  usually  repairing  the  use  of  computers.  Analytic  solutions  are  to  be 
preferred  if  feasible.  Sometimes  part  of  a  problem  can  be  solved  analytically 
end  then  the  intractable  parts  processed  by  simulation.  A  detailed  search  for 
shortcuts  and  approximations  before  proceeding  with  the  actual  simulation  will 
often  result  in  sizeable  savings  of  computer  time. 

4.  The  engineer  will  usually  be  forced  to  simulate  with  one-dimensional 
spectral  densities  because  there  are  so  few  reliably  measured  two-dimensional 


0^  values  is  shown  in  Fig.3» 
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Fig.  4. 


Cq  values  from  the  simulations  plotted  versus  Reynold's  number. 
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o 

Fig,  5  •  Cjj  from  simulations  plotted  on  normal  paper. 
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spectral  densities.  The  formulas  for  2-D  spectra  were  included  in  the  hope 
that  satisfactory  directional  data  will  soon  be  available. 

5.  The  linearization  of  the  drag  force,  or  its  approximation  by  poly¬ 
nomials,  is  particularly  useful  in  many  problems.  However,  in  some  applications, 
it  is  just  as  easy  to  work  with  the  v|v|  term  directly. 
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APPENDIX  I.  -  FOURIER  TRANSFORM  OF  A  WEIGHTED  DIRAC  COMB 
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where 


A«  =  |  <*„  +  *-„>  <100> 

B„  =  |  <*n  '  *.„>  <101> 


F  =  1/2  A  t 

Ihls  definition  results  in  the  symmetries 
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The  previous  equation  reduces  to 
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APPENDIX  II.  -  INPOT-OUTPUT  RELATIONS  FOR  MULTIPLE  SIMULATIONS 


Using  the  notation  of  eq.  (42)  and  following 
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One  summation  is  eliminated  by  the  independence  of  the  inputs.  The  summand  is 
zero  unless  j  =  j'  .  By  definition 
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If  (112)  is  substituted  into  (110)  and  it  is  assumed  the  functions  involved 
permit  interchange  of  integration  order,  then 
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(The  over  bar  indicates  complex  conjugation.) 

Again,  by  definition  of  Cmr(T),  analogously  to  (112) 


Cmr<T>  = 
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(114) 


If  (114)  is  compared  with  (113),  it  can  be  seen  that 


8mr<f>  =  Z  3xlXl(f) 

jsl 


(115) 


If  the  inputs  are  white  noise,  the  input  spectral  densities  are  unity  and 

■w<f>  =  Z  Vf)  W77  (11£> 

Jsl 

Since 


Sr<f>  =  »ra<f> 

(116)  determines  s_-(f)  for  all  m  and  r  between  1  and  M. 
mr 
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APPENDIX  III.  -  NOTATION 


The  following  symbols  are  used  in  this  paper: 

a  r  horizontal  component  of  acceleration; 

an  r  weighting  coefficient  in  a  digital  filter  (see  eq.(l?)); 

2 

a  z  increment  in  S(f)  (see  eq.(7)  ff.); 

A  z  a  constant  in  the  Bretschneider-Pierson  spectral  density  formula 
(see  eq. (10)); 

An  :  coefficients  in  Fourier  transform  of  a  digital  filter  (see  eq. (30); 
B  z  a  constant  in  the  Bretschneider-Pierson  spectral  density  formula 
(see  eq.(10)); 

Bn  z  coefficients  in  Fourier  transform  of  a  digital  filter; 
c  -  coefficient  in  optimal  linear  approximation  of  drag  force 
(see  eq.  (72)); 

cn  z  coefficients  in  optimal  polynomial  fitting  to  drag  force; 

Cq  r  drag  coefficient; 

Ctf  Z  inertial  coefficient; 
d  s  water  depth; 

D  =  pile  diameter; 
f  :  frequency  in  cycles  per  unit  time 

A 

fn  z  midpoint  of  the  n-th  increment  of  frequency  (see  eq.(2)); 

F  z  cutoff  frequency  -  the  highest  frequency  important  in  the  response 
or  a  frequency  such  that  s(f)  ssO  for  f  >  F; 
g  =  gravity; 

h  =  space  lag  in  x  direction  (see  eq.(55)); 

Aw  =  imaginary  part  of  z; 
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wave  number; 

kernal  in  a  linear  operator  (see  eq.(l8)); 

wave  number  in  radians  per  unit  length  corresponding  to  frequency 

A 

f  (see  eq.  (4)); 

Fourier  transform  of  k('i'),  K(f)  is  called  the  system  function; 
mean  velocity  in  eq.(83); 

number  of  increments  in  simulation  by  superposition; 

quantity  to  be  minimized  in  least-square  fitting  (see  eq.(?2)); 

coefficients  in  the  normal- to- gamma  transformation  (see  eq. (50); 

real  part  of  z; 

skewness; 

spectral  density; 

spectral  density  of  x(t); 

cross- spectral  density  of  x(t)  and  y(t); 

cross-spectral  density  between  ym(t)  and  yr(t)  in  the 

development  following  eq. (42); 

CumnulatLve  spectrum  (see  eq.(6)); 
time; 

horizontal  component  of  velocity; 

root-mean- square  velocity,  or  standard  deviation  of  the  velocity; 

specific  weigh v  of  water; 

one  of  the  horizontal  coordinates; 

v/v^  (see  e<**(79)  ff.); 

variate  in  standard  normal  formulas  (eqs.(84)  and  (85)); 
general  input  for  linear  operators; 

n-th  element  in  input  sequence  for  a  digital  filter  (see  eq.(l?)); 
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Xn  : 


y  = 


yk  = 


z  Z 
S(t)  = 

A^n  z 
At  = 
z 

>1(  t)  = 
VbM  = 
ftt)  = 

5  = 
9  Z 

K  = 

cr  - 
= 

IpU)  : 

$- 
$(x)  ; 


coefficients  in  the  gamma-to-narmal  transformation  (eq.(52)); 
one  of  the  horizontal  coordinates; 

k-th  element  in  the  output  sequence  of  a  digital  filter 
(see  eq. (1?)); 

coordinate  measured  vertically  upward  from  sea  floor; 

Dirac  delta  function; 

n-th  increment  of  frequency  (see  eq.(l)); 

time  increment  in  application  of  a  digital  filter; 

m-th  increment  of  6  ; 

sea  surface  elevation  above  still  water  level; 

A 

skewed  version  of  l^(t); 

simulated  sea  surface  elevation  above  still  water  level  at  time  t 

direction  variable  in  a  directional  spectrum; 

wave  phase  angle  in  eqs.(?0)  and  (71); 

midpoint  of  the  m-th  increment  of  B  ; 

standard  deviation  of  velocity  r  vma; 

wave  force; 

standard  normal  probability  density  (eq. (84) ) ; 
random  phase;  and 

distribution  function  for  a  standard  normal  variate  (eq.(85)). 


